#!/usr/bin/perl
# 2014/01/30 fomartの変更に対応。このプログラムで閾値をコントロールするようにする
# 2011/02/16 significantファイルを本当のbed formatに変換

use strict;
use warnings;
use IO::File;
use Getopt::Std;
use Carp qw(croak);
$| = 0;


### 0を指定するとfilter無し
if(@ARGV < 2 or @ARGV % 2 == 1 or $ARGV[0] eq '--help'){
	die "Usage : $0 -i [file] -L [fold change:default 2] -P [Pvalue: default 1e-10]\n";
}

my %opt;
getopts("i:L:P:", \%opt);
my $file = $opt{i};
my $Threshold_fold_change = 4;
if(exists $opt{L}){
	$Threshold_fold_change = $opt{L};
}
my $Threshold_Pvalue = 1e-40;
if(exists $opt{P}){
	$Threshold_Pvalue = $opt{P};
}

##chr	start	end	MAX_pos	treatment	control	SUM_treatment	SUM_control	fold-change	P-value	Q-value
#I	50	8441	311	846.83	733.3431708	22139	19618	1.12850443470282	2.57288042257575e-35	4.41032080854536e-37
#I	18173	18577	18342	586.7	511.2743763	904	823	1.09842041312272	0.0242215824293514	0.000194165157355866
#I	60604	60944	60806	509.62	284.2685034	701	407	1.72235872235872	2.42875798444658e-19	3.26088818148695e-21
#I	65846	66222	66042	787.09	311.347631	1025	446	2.29820627802691	2.43461898915709e-53	4.98291156387022e-55


my $COUNT_TOTAL = 0;
my $COUNT_final = 0;
my $COUNT_passFold = 0;
my $COUNT_passQval = 0;
my $fh = IO::File->new($file) or die "cannot open $file: $!";
while($_ = $fh->getline()){
	s/\r?\n//;
	if(m/^#/){
		next;
	}
	my $Flag_pass = 0;
	$COUNT_TOTAL++;
	my ($chr, $start, $end, $MAX_pos, $treat, $control, $SUM_treat, $SUM_control, $fold_change, $Pvalue, $Qvalue) = split /\t/;
	if($Threshold_fold_change == 0 or $fold_change > $Threshold_fold_change){
		$Flag_pass++;
		$COUNT_passFold++;
	}
	if($Threshold_Pvalue == 0 or $Pvalue < $Threshold_Pvalue){
		$Flag_pass++;
		$COUNT_passQval++;
	}
	if($Flag_pass == 2){
		$COUNT_final++;
		my $peakName = sprintf "peak%d; %.3f; %.3e; %.3e", $COUNT_final, $fold_change, $Pvalue, $Qvalue;
		print "$chr\t$start\t$end\t$peakName\t0\t+\n";
	}
}
$fh->close();


my $message =<<"EOF";
Total peak candidate : $COUNT_TOTAL
Passed fold-change filter (threshold $Threshold_fold_change) : $COUNT_passFold / $COUNT_TOTAL
Passed Pvalue filter (threshold $Threshold_Pvalue) : $COUNT_passQval / $COUNT_TOTAL
Final peaks : $COUNT_final
EOF

	warn "$message\n";